Nonequilibrium pairing instability in ultracold Fermi gases with population imbalance 
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We present detailed numerical and analytical investigations of the nonequilibrium dynamics of 
spin-polarized ultracold Fermi gases following a sudden switching-on of the atom-atom pairing cou- 
pling strength. Within a time-dependent mean-field approach we show that on increasing the imbal- 
ance it takes longer for pairing to develop, the period of the nonlinear oscillations lengthens, and the 
maximum value of the pairing amplitude decreases. As expected, dynamical pairing is suppressed 
by the increase of the imbalance. Eventually, for a critical value of the imbalance the nonlinear os- 
cillations do not even develop. Finally, we point out an interesting temperature-reentrant behavior 
of the exponent characterizing the initial instability. 

PACS numbers: 03.75.Ss, 03.75.Kk. 



I. INTRODUCTION 

One of the new exciting avenues that can be explored 
in the study of many-body properties of cold atomic 
gases [II 12] is the nonequilibrium dynamics following a 
sudden quench. Present-day technology allows to change 
the coupling constants on such short time scales that 
it is possible to explore the regime where the many-body 
system is still governed by a unitary evolution but with 
nonequilibrium initial conditions. Time-dependent cou- 
plings can be realized, for example, by varying the inten- 
sity of the laser that fixes the amplitude of an optical lat- 
tice or by changing the atomic scattering length through 
sweeping an external magnetic field across a Feshbach 
resonance. This problem, which has attracted a lot of 
attention recently |l[51ll[71iaiilini[IIJ[nj[I31[TllI51 
P^6, 17, 18 , is what we consider as well. 

Our work is inspired by Refs. |5] and |7] that deal with 
the study of the dynamical pairing instability in cold 
atomic gases after a sudden switch of the attractive in- 
teraction at times shorter than the quasiparticle energy 
relaxation time. Barankov et al. [5], starting from a nor- 
mal state, showed that after the quench the system is 
unstable. Pairing correlations initially build up expo- 
nentially in time and then oscillate taking the form of 
soliton trains. If the system before the quench is in an 
equilibrium BCS state, and the quench is performed by 
changing abruptly the pairing coupling, then the station- 
ary state can show a constant (but reduced) gap or can be 
gapless [m [T9j . A classification of the allowed nonequi- 
librium behaviors arising from different initial conditions 
has been presented in Ref. 20, To date there are no ex- 
periments on the non-adiabatic switching of pairing in 
fermion condensates. A proposal to detect signatures 
of nonequilibrium dynamics using radio-frequency spec- 
troscopy has been put forward recently [16 . 

Along the lines of these previous works (see also 
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Ref. [21]), in the present paper we study the pairing in- 
stability in a two-component ultracold Fermi gas with 
unbalanced spin populations after a sudden switch of the 
attractive interaction between the two fermion species. 
As is well known since the early days of superconductiv- 
ity [221 [231 [21] , an imbalance in the number densities of 
the two species tends to suppress pairing. Unbalanced 
Fermi gases [25] are currently attracting a great deal of 
experimental and theoretical interest. One of the aims is 
to detect exotic paired states [211 [2S1 1271 12H] that have 
been elusive so far in conventional solid-state systems. 
Fermi gases with population imbalance have been real- 
ized in a series of experiments [2S1 [SHI [33 132] • The equi- 
librium phase diagram has been worked out in great de- 
tail (see, for example, Refs. [331 [31 [33 [33 [331311 and 
references therein) and a very rich scenario has emerged. 
However, despite the tremendous effort that has been 
devoted to understand equilibrium phases, nothing is 
known yet about the out-of-equilibrium properties of 
these system. Here we address this question for the first 
time. As a first step we analyze the instability of a normal 
partially spin-polarized Fermi gas with respect to s-wave 
pairing which leads to nontrivial results. Guided by the 
body of knowledge acquired in the study of the equilib- 
rium case, one can look also for instabilities towards more 
complex paired states that we leave for future study. 

The time scales that are relevant to the present prob- 
lem |12j are the quasiparticle Landau Fermi-liquid life- 
time Tel, the time over which the oscillations of the 
pairing function develop and evolve [5], and the charac- 
teristic time To over which the coupling is switched on. 
We are interested in the regime when the inequalities 
To < TA < t < Tci hold. 

The paper is organized as follows. In the next Section 
we first introduce the model Hamiltonian that we use to 
describe the system of interest. In Section [HA] we discuss 
the mean-field decoupling used to study the time evolu- 



tion, while in Sect. II B we carefully describe the initial 
state to which the quench is applied. The resulting equa- 
tions can be analyzed both numerically and analytically. 
In Sect. |IIl"S] we present our numerical simulations of the 
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time-dependent mean-field equations and discuss their 
main features. In Sects. IIIB and IIIC we present some 



commuting terms Ti.MF(t) 



analytical results for the short- and long-time properties 
of the quantum evolution. In Sect.lTVlwe summarize our 
main conclusions. Finally, Appendix|A|contains more de- 
tails on the numerical simulations of the time-dependent 
mean-field equations, while Appendices [B| and [C] contain 
some details of the calculations presented in Sect. |IIIB[ 



II. THE MODEL 



The time-dependent BCS Hamiltonian is defined as 



k.a 



k.k' 



(1) 

In this equation cj,^ (cka-) creates (annihilates) a fermion 
with momentum k {h = 1) and spin cr =1,1 (hyperfine 
state label) . The number of particles with spin a is 
fixed during the time evolution and thus we do not need 
to introduce chemical potentials for each spin species [55] . 
Given N^, the equilibrium Fermi energies epj and epj. of 
the noninteracting system at zero temperature are fixed. 
The summations in Eq. (jlj are carried out over a shell 
of energies of thickness 2wd around the Fermi energies, 
where wd is an effective ultraviolet cutoff frequency |40j . 
We assume that the Fermi energy mismatch, 6fi = ept ~ 
epj, is smaller than wq. For convenience we measure 
all the energies from epj, and approximate the parabolic 
dispersion Sk with a sequence of iV ^ 1 equally spaced 
levels £fe in the range [— wd , '^d] i where fc = 1 . . . A'^ is a 
scalar label. The level spacing is Se — 2uj-[)/{N — 1) and 
the density of states is l/6e. 

The coupling g{t) is zero if i < and is switched on 
to a constant negative value —g during a time interval 
< t < tq. Since we focus on the non-adiabatic evolu- 
tion (to <^ ta), we approximate g{t) « ~'gQ{t), where 
0(x) is the Heaviside step function. It is worth to notice 
that if the switching on of the interaction is too fast, the 
gas becomes overheated and the time-dependent coupling 
induces two-particle scattering. However, as discussed in 
Ref. [121 a time window for tq exists in which the con- 
straint for avoiding the overheating is compatible with 
that of a sudden switching-on of the interaction. 



A. Time-dependent mean-field tlieory 

As discussed in Ref. I12| the nonequilibrium evolution 
of the fermion system can be analyzed within a time- 
dependent mean-field theory. To this end we introduce 
the pairing function A(t) = g^^ic-kiCk^), where the 
average is taken over the quantum state of the system at 
time t. After the mean-field decoupling is performed, the 
BCS Hamiltonian ([T]) reduces to a sum of time-dependent 



Efc^MpW' where 



(7 

... (2) 

Within the mean-field approximation the Hilbert space 
to study the time evolution of the system is the tensor 
product of N Fock spaces with at most two particles in- 
stead of the larger Fock space with at most 2N particles. 
There are only four states in the two-particle Fock space 
built with the single-particle orbitals: the vacuum state 
|0), a fully-occupied state |2) with two particles, and two 
singly-occupied states |t) and ||) labeled by the spin of 
each unpaired fermion. Writing the Fock basis in this 
order, the matrix Ti-Mpit) within a block with a given k 
is 






-A*(t) 








A(t) 
















£k 

















(3) 



The Hamiltonian decomposes into four blocks along the 
diagonal. The last two blocks are one-dimensional and 
determine the free evolution of the unpaired states, as 
these states cannot be coupled to the |0) ® |2) conden- 
sate sector due to the Pauli- blocking effect. The two- 
dimensional block represents a Cooper pair, where the 
vacuum |0) is coherently coupled to the doubly-occupied 
state |2). The coupling is due to the pairing term c[.|clj.| 
that does not conserve the number of particles within the 
subspace. 

Since it is important to include the case where the 
fermions can be excited out of the condensate into un- 
paired states by incoherent thermal processes, a wave 
function is not appropriate to treat the evolution of the 
two-particle system. To treat this problem we use a sta- 
tistical matrix defined as 

p^''\t) = (i - - p^ki)[uk{m + vkit)\2)] 
X mulit) + {2\ii{t)] +miT)(Ti +P-fciii)ai . 

(4) 

The probabilities pk] and p^ki take into account the ther- 
mal excitation of particles out of the condensate. We 
remark that each pure state that enters the construc- 
tion of the statistical matrix has to be normalized, i.e. 
\uk{t)\^ + \vk{t)\' = l. 

Both the Hamiltonian ([3| and the statistical matrix Q 
are block-diagonal and the condensate sector evolves in- 
dependently of the other states, according to idtp^^\t) = 
[^mf(^)' P^'^H^)]- We can define an effective Hamilto- 
nian restricted to the condensate sector and an ef- 
fective state vector |Sfe(t)) = Uk{t)\Q) + Vk{t)\2), with 
Uk{t) = Uk{t)(\ -Pk-\ -p-fex)^/^ and Vk{t) = Vk{t){\ - 
Pk\ — p-ki)^^"^ ■ The statistical matrix projected onto 
the condensate sector then reads pi''\t) = |Sfc(t)) (Sfc(t) |. 
The pure-state form of the projected statistical matrix is 
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preserved by the time evolution. This implies that the 
effective, non-normalized state vector belonging 
to the condensate sector |0) ® |2) is sufficient to describe 
the time evolution. 

The state vector |Sfe(i)) evolves according to the norm- 
preserving effective Schrodinger equation idt\^k{t)) — 
'Hi''\t)\'E.k{t)) , and the coefficients Uk{t) and Vk{t) obey 
the time-dependent Bogolubov-de Gennes equations 
(BdGE) 



idt 



Vk{t) 
Uk{t) 



Sk -A{t) 
-A*(<) ~ek 



Vk{t) 
Ukit) 



(5) 



The total Fock space for (at most) 2N particles is then 
defined to be the tensor product of the two-particle 
spaces and the statistical matrix is p = (^^.p^'^-'. If 
an operator Ok has support within the condensate sec- 
tor of the k space, its expectation value Tr[/3C'fe] can be 
computed using the effective state vector only and reads 
{Ek{t)\dk\Sk{t)). The BdGE have to be solved together 
with the self-consistency condition 



A{t) ^ gJ2uUt)vk{t) . 



(6) 



variable of fc. As a consequence we can take as initial 
conditions 



Uk = ywvyw^ 

Vk = exp [i(l)k] ^/fk^fki ■ 



(9) 



A non-zero temperature or a finite value of the imbal- 
ance are sufficient to produce a non-zero initial pairing 
amplitude |A(t = 0)|, which is very small because of the 
randomness of the initial phases. 



III. RESULTS 

In this Section we discuss our results for the time 
dependence of the pairing A(t) and the distribution of 
paired particles nk(t) as functions of spin imbalance, tem- 
perature, and initial conditions. The numerical results, 
obtained through integration of the BdGE, will be sup- 
plemented by analytical results obtained in the short- 
time and stationary regimes. 



Numerical solution of the BdGE 



B. The initial state 

The BdGE in ^ must be accompanied by some initial 
conditions Uk — Uk{t = 0) and Vk = Vk{t = 0). The 
initial conditions thus describe the state of the system 
just before the quench is applied at time i — > O"*" . We have 
chosen initial conditions corresponding to the equilibrium 
configuration of the Hamiltonian ([2| at temperature 9 
and g = 0. We compute the partition function Zk of the 
A;-th subsystem in the grand-canonical ensemble 

Zk — l-|-e^'^^^^''~^T-A'l) _j_g-/5(£fc-A'T) -)_g-/3(efc-AH) ^ (7) 

where (3=1/9 (fee = 1), ^^'^ chemical 
potentials for the two spin species, and the difference 
/i| — fii is equal to the Fermi energy mismatch dp,. 
The probability to find the system in the state |2) is 

\Vk\' = (2|pW|2) = ^exp[-/3(2£,-^T-/ii)] 

= fk]fkl , (8) 

with = {1 + exp[/3(£fc - Sp)]}^^ and fki = [1 + 
exp(/3efc)]~^. Similarly, the probability to find the sys- 
tem in the state |0) is \Uk\^ = (1 - /fcT)(l - /fcj)- The 
probability | C/fc P -t- 1 Vfc p to find the fc-th subsystem in the 
condensate sector is smaller than unity: because of ther- 
mal excitations there is a finite probability that the k-th 
subsystem is occupied by an unpaired fermion. It is easy 
to see that the expression |ufe(t)p -I- |wfe(t)p is constant 
in time. 

Since at times t < the system is noninteracting, the 
phase (pk of the coherence (2|p''^^|0) = U^Vk is a random 



We now turn to the presentation of the numerical so- 
lution of the BdGE ([5| with initial conditions given in 
Eq. ^ . In what follows we use as unit of energy the real 
quantity Aq defined by the solution of the equilibrium 
BCS self-consistency equation g^ki^k + Aq)~-^/^=2. 
This choice of the energy scale then fixes the value of 
g. Frequency and time scales are defined accordingly. To 
solve the BdGE we have used a fourth-order adaptive- 
stepsize Runge-Kutta algorithm, with a maximum rela- 
tive error of 10~^ per time step. A typical time step is 
10~^ — 10~^, but a smaller time step of order 10~^ is used 
near the initial instability of the BdGE (see below). The 
integration of the BdGE up to t,„ax = 300 takes less than 
10 sees on a desk PC. 

In Fig. [l] we show some representative results of the 
solution of the BdGE for iV = 10^, wd = 5.0 and g ~ 
4 X 10^'^. We choose three initial states with different 
imbalance Sp at a temperature 9 = 10^^. Each proffie is 
obtained with a random realization of the initial phases 
(j)k that we take as uniformly distributed in the interval 
[0,27r]. 

Three time regimes are evident for each value of the 
imbalance Sp in Fig. [l] (i) a very short initial transient 
[0, iin] where the pairing amplitude increases by several 
orders of magnitude, as will be clarified in Sect. IIIB (ii) 
a time interval [tin, t] in which the growth of |A(t)| is ex- 
ponential in time, |A(t)| = 77 exp (7^) (r will be hereafter 
referred to as "time lag" , following the jargon introduced 
in Ref. [5|); and (iii) a time interval where undamped, 
nonlinear oscillations of |A(t)| occur. 

Several observations are in order at this point. On 
increasing the imbalance Sp the exponent 7 of the expo- 
nential growth in region (ii) decreases and the time lag 
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FIG. 1; (Color online) Modulus jA(f)j of the pairing function 
(in units of Aq) as a function of time (in units of I/Aq), ob- 
tained by solving the BdGE. From top to bottom the value of 
the initial Fermi-energies mismatch increases as Sjj, — 0.0, 0.5, 
and 0.75. The left panels show a zoom of the initial linear 
instability region in the range t < t, t being the time at 
which |A(t)| has its first peak. The time interval [0,tin] is the 
transient discussed in (i) in Sect. Ill A The thick dashed lines 
are linear fits in the range 0.3r < t < 0.8t: the slope of each 
dashed line gives 7, while the extrapolation to t = gives 
77. A computer precision of 10~^^ is reached for S/i = 0.75 
(bottom panel) and t < 10.0, where fluctuations due to the 
numerics begin to appear. 



r increases (i.e. it takes longer for pairing to develop), 
the period of the nonhnear oscillations lengthens, and the 
maximum value of the pairing amplitude decreases. As 
expected, dynamical pairing is suppressed by the increase 



of the imbalance. In Sect. IIIB we prove that dynamical 
pairing is wholly suppressed at a critical value SfXc of the 
imbalance. It is hard to verify this assertion numerically 
because at large imbalance the initial pairing |A(t = 0)| 
becomes comparable to the computer accuracy. 

To test the robustness of the profiles shown in Fig. [T] 
against changes in the initial conditions we have solved 
the BdGE with several different choices of the initial ran- 
dom phases. The results of this statistical analysis are 
reported in Appendix [Xj where we show that the am- 
plitude of the pairing is essentially independent of the 
particular realization of the random phases. 

In Fig. [2] we show the distribution of condensed parti- 
cles 
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FIG. 2: A three-dimensional plot of the quantity 7ifc(f) — nfe(O) 
as a function of energy Ek and of time t. In the top-right panel 
we show Ukit*) — nk{0) as a function of Ek at a time instant 
t* where the pairing amplitude is maximal. In this figure 
= 0.5, as in the central panel in Fig. [l] 



measured from its initial value nk{0), as a function of 
energy Sk and time t. As a function of time, the quan- 
tity nk{t) is always nearly equal to its initial value nfe(O) 
except in close proximity to the maxima of the pairing 
amphtude |A(t)|. As time evolves, nk{t) — nk{0) pulses in 
synchronism with the nonlinear oscillations of the pair- 
ing function. Close to a time t* at which the pairing 
amplitude is maximal, nk{t) — nk{0) exhibits a peculiar 
structure (see top-right panel in Fig. |2]) . We in fact see a 
downward peak in the region below the Fermi surface of 
the minority-spin component and an upward peak, equal 
in size to the downward one, located above the Fermi 
surface of the majority-spin component. In between the 
two peaks we recognize a region of extension where 
pairing is suppressed since the condensate sectors |0)® |2) 
are almost entirely depleted, i.e. \uk{t*)\'^ + \vkit*)\'^ ~ 
for Sk G [0,(5^]. The two peaks indicate that particles in 
the condensate are transferred across the Fermi surfaces 
of the two populations. This phenomenon is reminiscent 
of what happens in conventional BCS equilibrium super- 
conductivity. 

In what follows, we show that the different regimes 
of the initial onset of the pairing instability and of the 
nonlinear oscillations are amenable to an analytical treat- 
ment. In particular, in Sect . [HTB] we solve by means of a 
linear-stability analysis the time-dependent BdGE in the 
time interval [0, t] (regions (i) and (ii) introduced above). 
In Sect. |III C| we discuss the stationary limit within the 
general theoretical framework that was earlier developed 
in Refs. lHj and [20, for the unpolarized case. The main 
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result of these two sections is a complete analytical pre- 
diction of the solutions of the time-dependent BdGE. 



B. Analysis of linear instability 

The initial build-up of the pairing instability can be 
studied by means of a linear-stability analysis, along the 
lines of what was earlier done in Ref . 5 for the unpolarized 
case. 

It is convenient to introduce the following definitions, 
corresponding to a free evolution of each Cooper pair, 



times t > tj, 



Uk{t) 
Vk{t) 

m 



Uk 

k 



(11) 



where Uk and Vk are the initial values in Eq. ([9| . Without 
loss of generality, we can write any solution of the BdGE 
in the form Uk(t) = Uk{t) + 5uk{t) and Vk{t) = Vkit) + 
Svk{t). We choose Suk{0) — 6vk{0) = so that the initial 
conditions are still given by Wfc(O) = Uk and Vk{0) = Vk- 
Inserting these definitions into the BdGE we obtain the 
equations of motion for the corrections 5uk{t) and 5vk{t), 



idt 



Suk{t) 
6vk{t) 



-EkSukit) - A*{t)[vk{t) + Svk{t)] 
-A{t)[uk{t) + Suk{t)] + ekSvkit) 



^ (12) 
We solve Eq. (12 1 in a time interval < ^ ^ ^ defined 



by the hypotheses 



(i) |A(t)|»|A(i)| 

(ii) \Svk{t)\ <t: \vk{t)\ 



(13) 



These hypotheses mean that after an "instability time" 
iin the pairing function built up by the corrections 5uk 
and Svk is much larger than the pairing due to the un- 
perturbed functions Uk and Vk- The first hypothesis 
is fulfilled if the initial state is weakly paired, i.e. if 
A(0) = g\ J2k ^k^kl < Aq. The second hypothesis guar- 
antees that the corrections are much smaller than the 
unperturbed functions, so that we can neglect the nonlin- 
ear terms 6ul.Svk in the pairing function. The nonlinear 
terms become important only after a "nonlinearity time" 
<r. 

The time evolution in the interval [ti„ , r] is ruled by 
the linear ordinary differential equation 



idt 



5uk{t) \ _ 
Svk{t) 



-EkSukit) - SA*{t)vk{t) 
--SA{t)uk{t) + EkSvkit) 



(14) 



where 6A{t) = gJ2kKit)SMt) + Sul{t)vk{t)]. This 
equation does not allow us to trace the nonlinear evo- 
lution in the interval [0,tin]. We only need to assume 
that 5uk(tin), 5vk(tin) and SAiti^) are non-zero and we 
write the following Ansatz for the solution of Eq. (14) at 



dA{t) 
6uk{t) 
Svk{t) 



-iC(t-tin) 



5A{tin) 



(15) 



Here we have introduced a complex instability exponent 
( = io + i^. Inserting the Ansatz ( 15 1 in Eq. ( 14 ) one can 



easily obtain [12] the following "consistency relation" for 
the instability exponent ^, 



E 



\Uk\'-\Vk\' 



1 

9 







(16) 



This equation is identical in form to Eq. (18) in Ref. [T^l 
but here the solution ( = C{S^, 9) depends on two physi- 
cal parameters: the imbalance (5/i and the temperature 9 
(rather than only on temperature, as in the unpolarized 
case). For J/^ = we recover the results in Fig. 10 of 
Ref.m 

In Fig. [3]we show the imaginary part of the solution of 
Eq. ([T6| in the (5/^, 9) plane. To solve Eq. ( 16 ) we have 
minimized the square of the l.h.s. with respect to the 
two parameters uj and 7. The minimum of the square is 
just the value where the l.h.s. vanishes. Several observa- 
tions need to be done on Fig. [3] To begin with, there is 
a critical line in the ((5/i, 9) plane above which no insta- 
bility develops, i.e. 7 = 0. The imaginary part 7 of the 
instability exponent decreases monotonically as a func- 
tion of On the contrary, 7 depends monotonically 
on temperature only \i Sji < S^r — 0.7. In this case 7 
decreases if 9 increases, while the opposite behavior hap- 
pens if Sfi > Sfir and the temperature is low. The latter 
region of the (5/i — 9 plane appears as a re-entrance in 
the bottom panel of Fig. [3] In this region an increase 
in temperature allows the system to sustain pairing even 
in the presence of a larger maximum imbalance. This is 
reminiscent of a similar re-entrant behavior obtained in 
the equilibrium case by Sarma [23. . In that case, how- 
ever, the author found the existence of a more stable 
phase characterized by the absence of re-entrance. The 
calculations in Ref. [53] are equilibrium calculations per- 
formed within a grand-canonical ensemble and thus do 
not rule out the possibility of a re-entrance in the "phase 
diagram" of Fig. |3] for the out-of-equilibrium dynamics. 

In some limiting cases it is possible to extract analyt- 
ically the solution of Eq. (16). In the thermodynamic 
limit, defined by letting N —^ 00 while keeping Aq and 
t^D fixed, Eq. ( 16 1 reduces to jlT] 



de 



2e-uj 



1 



[l-/T(e)-/i(£)] 



5e 







de j-^ , [1 - h{e) - /x(e)] = , (17) 

(2£ — ti>)^ + 7-^ 



where the real and the imaginary part have been written 
separately. The Fermi functions /a-(e) weigh the states 
that take part in the pairing process. The states in which 
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FIG. 3: Top panel: a plot of the imaginary part 7 = "/{Sfi, 9) 
of the instability exponent as a function of Sjj, and 6. Bottom 
panel: contour plots corresponding to the top panel. The 
thick solid line shows the points of the {Sfi, 9) plane where 
the 3D profile in the top panel intersects the 7 = plane 
(in actuality this curve has been calculated for 7 = 10~^ for 
numerical reasons). The re-entrance described in the main 
body of the text is clearly visible. 



value is consistent with the Thouless criterion for super- 
conductivity [42] )• We remind the reader that supercon- 
ductivity is suppressed by the apphcation of a Zeeman 
field larger than the critical Clogston-Chandrasekhar 
value [3^, which translates into a critical imbalance 
^^J'CC — %/2. Note also that the transition (18 1 from 
the paired to the unpaired regime is continuous with a 
singularity in the derivative, as in a phase transition of 
the second kind. Subleading terms in the asymptotic ex- 
pansion in powers of 1/wd are presented in Appendix [B] 
and do not modify the key features of Eq. (18 1. 

We now study Eqs. ( 17l for small but finite in order 
to determine the value of the imbalance 5^^ above which 
the dependence of 7 on ceases to be monotonic, i.e. 



7(6*, (5^) > 7(0,5^) if Sfi>Snr 



(19) 



We expand 7 near ^ = 0, j{6, Sfi) — ■y{9 = 0, S^) + 
6ji{6fi) + 6'^72((5/i) -I- 0{6^). A similar expansion is writ- 
ten for Lu(9,S^). The integrals involving the Fermi func- 
tions in Eqs. (17) can easily be computed up to second 
order in 9 using the Sommerfeld method, as briefly out- 
lined in Appendix [Cj In the limit uj-d 3> 1 we obtain 
7i(i5/i) — and 



(20) 



We see that 72 > for 6fi > %/2/2, i.e. S^x^ = V2/2 
and 7 increases quadratically with temperature. In Ap- 
pendix [C] we report an expression for Sji^ that is correct 
up to second order in I/uj-q. 

Before concluding this section, we would like to men- 
tion that the existence of a re-entrance for > 1, i.e. 
d'^5fi/dd'^\.y=o > 0, can be proven by arguments similar 
to those that led to Eq. (l20|. 



C. Analysis of the pairing oscillations 



there is a high probability to find an unpaired electron 
are effectively removed from the system. This is most 
clearly seen at 6 = 0, where the Fermi functions become 
sharp steps and 1 — /|(e) — /^(e) = 6(e — (5/i) — 6(— e), 
thus excluding the interval [0,(5^] from the integrations 
in Eq. (17). We see that the exclusion of some fermions 



from the pairing must lead to a decrease in the expo- 
nent 7 of the instability, or equivalently in the maximum 
amplitude A_|_ of the oscillations. In the zero tempera- 
ture 9 = case (see Appendix [b]) , after performing an 
asymptotic expansion in powers of 1 /ujd we find that the 
solution of Eqs. (17) is 



j{9 = 0,6fi) = ^l-Sfi^ , 



(18) 



for < Sji < 1. We thus see how an imbalance larger 
than (5/ic = 1 inhibits the development of pairing (this 



In this section we focus on the oscillatory dinamics of 
the pairing function, shown in the right panels of Fig. [T] 
We follow Refs. [TTJ [201 and [321 and use the formalism of 
the so-called Lax vector that allows an implicit analytical 
solution of the BdGE. 

The Lax vector L{w) is a three-dimensional vector 
whose components are rational polynomials of an aux- 
iliary complex variable w and is defined as |43j 



L{w) 



E 



(21) 



Here z is the unit vector in the z-direction and Sk = 
(S^, Si, S^) is a. three-dimensional real vector whose com- 
ponents are defined by 5"^ - iS^ = U^Vk and 23^ = 
|Vfcp — |C/fcp. According to Ref. [TT|the asymptotic time 
evolution of the solutions of the BdGE can be predicted 
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by looking at the roots of In the limit N ^ oo 

almost all the roots of cluster together on the real 

axis. Few isolated roots with non-zero imaginary part 
define the frequencies that appear in the oscillations of 
A{t). 

The vectors {Sk,k — 1...N} can be interpreted as 
Anderson classical pseudospins [44 . Each fc-pseudospin 
represents the state of a Cooper pair and the initial 
state (C/fe, Vfc) can be formally mapped onto a pseudospin 
chain. In the case of the initial state written in Eq. ([9]), it 
is easy to see that a substantial probability to find a 
Cooper pair in the doubly-occupied state |2) corresponds 
to a very small probability | C/fc P to find it in the vacuum 
state 10). To simplify the expression of the Lax vector in 
Eq. (21 1 we introduce, however, a more stringent condi- 
tion. We take gN\UlVk\ ^ 1, i.e. we assume that the 
initial pseudospins are almost entirely aligned in the z 
direction. The Lax vector then becomes 



L{w) 



1 

9 



^ 2e, - 2w 



(22) 



For £k < 0, the pseudospins are aligned along the +z 
direction and represent doubly occupied states, while for 
< Ek < Sfi the norm of the pseudospins \Sk\ is negligi- 
ble and vanishes at zero temperature, and for Sk > Sfi the 
pseudospins are aligned along the —z direction and rep- 
resent vacuum states. The length of the fc-th pseudospin 
\Sk\ gives the probability that the fc-th subsystem is in 
the condensate sector |0) (E) |2). So the states that con- 
tain unpaired electrons correspond to pseudospins with 
smaller length. 

In our case it is easy to see that all the roots of \L{w) p 
in Eq. ( 22 1 are doubly-degenerate and are given by the 
solutions C of the consistency equation (16 1 and their 
complex-conjugates. At this point we remind the reader 
that in Sect. |IIIB we have found a single solution of 
Eq. ( 16 1 (illustrated in the top panel of Fig. [S]) with non 
zero imaginary part. This implies that the root diagram 
of [^(w)!^ in the complex plane contains two degenerate 
vertical cuts. 

The corresponding solution of the BdGE has the 
form PU| 



A(t) = A+dn((A+(t-to),fc) 



(23) 



with fc^ = 1 — A^/A'^. Here dn(a;, k) is a Jacobi elliptic 
function and the maximum amplitude of the oscillations 
A+ is equal to the imaginary part of the root of \L{w)\'^, 
which we have just shown to be equal to 7 = 9m C. 
The period of the nonlinear oscillations can be written 
in terms of the complete elliptic integral of the first kind 
K{x) as 



2 

AT 



(24) 



The parameter A_ is not fixed by this analysis and de- 
pends on the values of . The distribution of S'^ de- 
pends on the particular realization of the random phases 
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FIG. 4 

described in Sect 



A comparison between the results of the simulations 
(circles) and the analytical results 



III A 



of Sect. IIIB (linesJT All the numerical results shown are 
average values over fifty simulation runs, for e = 10^^ (see 
Appendix [A} Fig. |6|. Panel (a): the average maxima of A+ 
(circles) and the theoretical prediction given in Eq. | |18[ ) (solid 
line). Panel (b): the imaginary part 7 of the instability expo- 
nent (calculated as explained in Fig. [T| is shown to coincide 
with A+, the solid line being the A+ = 7 bisector. Panel (c): 
the amplitude A+ of the oscillations (circles) is shown as a 
function of the period T. The dashed (solid) line i s th e period 
T for A- = 10"^ (A- = 6 X lO"*), as from Eq. E^. 



so that we expect fluctuations in the value of A_ and 



T. 



In Fig. |4] we show that the numerical solutions of the 
BdGE illustrated in Fig. [l] agree very well both with the 
linear-instability analysis (Sect. IIIB I and with the anal- 
ysis based on the Lax polynomial. 



IV. CONCLUSIONS 

The presence of a population imbalance modifies dra- 
matically the dynamical pairing instability in a two- 
component ultracold Fermi gas when an atom-atom at- 
traction is suddenly switched on. In this work we have 
considered the case when the instability occurs via the 
s-wave pairing channel. We find that the dynamical in- 
stability is suppressed if the initial imbalance exceeds 
a critical temperature-dependent value, in analogy with 
what happens in the equilibrium situation. The exponent 
characterizing the linear-instability regime does not de- 
pend monotonically on temperature and shows an inter- 
esting re-entrant behavior in the temperature-imbalance 
plane. A similar behavior has been observed in equi- 
librium calculations since the early work of Sarma [23], 
though in that case the re-entrant behavior corresponds 
to a metastable state. In the dynamical situation the 
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variational principle on the grand-canonical thermody- 
namic potential is of course not present and such re- 
entrant behavior can indeed be observed. It is very in- 
teresting to understand how our findings show up in a 
radio- frequency spectroscopy measurement J6j . Another 
important aspect, which is currently under investigation, 
is to understand whether it is possible to access more 
exotic pairing states after a quench. 
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APPENDIX A: QUALITATIVE ANALYSIS OF 
THE NONLINEAR OSCILLATIONS 

The very regular shape of the nonlinear oscillations al- 
lows us to define an average period T and a maximum 
amplitude A+ for each simulated profile |A(<)|. In prac- 
tice, these quantities are calculated as follows. For each 
single realization of the random phases (/i/c we find the co- 
ordinates {^i, A_|_i}^"'"' of the first AAmax peaks by means 
of a cubic interpolation. Then we compute the averages 
T - E.(^.+i - t^)/(A/'max - 1) and A+ = E, A+i/A^^^x, 
and their standard deviations 5T and 5A+. 

In order to illustrate the robustness of the nonlinear 
oscillations shown in Fig. [T] we report in Fig. [5] an anal- 
ysis of their shapes and periods, as found for a total of 
thirty realizations of the random phases. We notice that 
the spread of both T and A-|_ diminishes with increasing 
imbalance, becoming comparable to the typical 5T and 
(5A4- that one finds in a single realization. That is, with 
increasing imbalance the quantities T and A-|_ become 
less and less dependent on the initial random phases. 

In Fig.[6]we present a more quantitative account of the 
effect of the random initial conditions on the magnitude 
of the fluctuations. We have computed the average T of 
the period and the corresponding standard deviation 5T 
over fifty realizations. We see that the relative fiuctu- 
ations 5T/T drop by one order of magnitude when 5ji 
spans the range [10"'^, 1] (see the top panel in Fig. pi. 
The average ST jT of the relative fiuctuations of the pe- 
riod increases instead by two orders of magnitude when 
spans the range [10^^,10^^], while it becomes com- 
parable to ^T/T for 8 [I > 10"^ 

Finally, in the bottom panel of Fig. [6] we illustrate the 
behavior of the relative fluctuations 5A+/A+ of the am- 
plitudes, which drop by three orders of magnitude when 
(5 /J, spans the range [10~^,1]. The average (5A_|_/A+ of 
the relative fluctuations of the amplitude remains of the 
same order of magnitude as (5A-|_/A+. 



APPENDIX B: CRITICAL IMBALANCE AT 
ZERO TEMPERATURE 

In this Appendix we determine analytically the critical 
imbalance b[i^ at zero temperature, defined by 7(0 = 
0, = 0. The value 10^ of lo at criticality has also to 



be determined to solve consistently Eq. ( 17 1 
Equation pT| at 6* = reads 



[(2c.D-cc>)2+72][(2a>D+c.)2 + 72] = 
= + yrT74)4[(2,5/i - u^f 7'] + 7' 



and 



(Bla) 



arctan[(2tjD — ~ arctan[(2(5/i — = 

= — arctan(cj/7) + arctan[(2a;D — ■ (Bib) 



We assume that 2(5 /ic > o^c, as is suggested by the numcr 
ical solution and also by the zeroth-order solution (18 1 
Then in Eq. (Bla I we put 7 = and obtain 



2(5/i 



■ LO 



In Eq. (Bib I we perform the limit 7^0 and find 
1 



(B2) 



1 



1 



1 



2a;D 



25/i 



2wd 



(B3) 



In deriving this result we have used that arctan (a/7) 
7r/2 — 7/a. 



Substituting Eq. (B2| into Eq. (B3) we obtain 

2 



^1 = 



1 + Jl + l/u;2 



Using this result back into Eq. (B2) we find 

2wd 



(B4) 



(B5) 



A second-order expansion of Eqs. (B4| and (B5l in pow- 
ers of 1/wd finally gives 



1 1 



1 LOy 



1 



3 1 

8^ 



(B6) 



In our computations — 5.0, so these second order 
corrections are of order 10^^ (cjc — 0.995 and (5/ic — 
0.985). 

Now we show that the slope of the curve 7(0 — 0, (5/i) 
is singular at the critical imbalance (5/ic and we find an 
asymptotic form for the profile. We make the Ansatz 
7 = a\/S^c ~ (5/^ and tu = ujc + k{S^c — S^j,). Substituting 
this into Eqs. (Bll and discarding powers of Sfi^ — 6fi 



higher than one we find that the Ansatz is consistent 
provided that 



= 2w, 1 



2 1 



1 1 



(B7) 
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FIG. 5: Qualitative analysis of fluctuations in nonlinear oscillations similar to those shown in Fig. [T] Imbalance increases from 
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APPENDIX C: SUBLEADING CORRECTIONS 

TO &yL, 

In the main body of the paper, immediately above 
Eq. ( 20 1 , we introduced an expansion of 7 in powers of 
temperature near = 0. The coefficient 71 (5/^) of the hn- 
ear term is identically zero, while the coefficient of 
the quadratic term has been given only for 00. The 

equation 72((5/ii) — defines the imbalance Sfi,- above 
which the dependence of 7 on ceases to be monotonic. 
In this Appendix we find the second-order corrections to 
the quantity dfj,r in powers of I/ciJd. 



To this end, we note that Eq. (17 1 can be written in 
the general form 



de9{e)[l~f{e)-f{e-^i)] = K 



(CI) 



with f{x) = l/{e^^ + !)• To compute 72 (^m) we need to 
expand this equation in powers of the temperature 9. In 
order to do so we follow a familiar Sommerfeld procedure: 



we perform an integration by parts in Eq. (CI I, expand- 
ing the primitive G{e) of g^e) in powers of 0. The Som- 
merfeld expansion of the integrals involving the Fermi- 
Dirac functions to order 0^ gives 



G(^d) 



G(-c^d) - G(0) 
'92G(e) 



e=0 



G(<5/i) 



e—Sfi 



K . 

(C2) 



For the first of the two Eqs. (17 1 the function G is given 
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by 



G{x) = - In [{2x - ujf + -/^ 



while for the second it is given by 

G(x) = — arctan ( — 

27 V 7 



(C3) 



(C4) 



perature 9 through the functions uj = Lo^O^SfXr) and 



7 = {6, 5 fly). We expand Eq. (C2| order by order in 



powers of 9 and subsequently in powers of I/wd- By 
imposing that 72(<5/ir) = we finally obtain 



V2 



1 1 

4^ 



(C5) 



We remark that G depends parametrically on the tem- 
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